The Local Variational Principle 
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A generalization of the Gibbs-Bogoliubov-Feynman inequality for spinless particles is proven and 
then illustrated for the simple model of a symmetric double-well quartic potential. The method 
gives a pointwise lower bound for the finite-temperature density matrix and it can be systematically 
improved by the Trotter composition rule. It is also shown to produce groundstate energies better 
than the ones given by the Rayleigh-Ritz principle as applied to the groundstate eigenfunctions of 
the reference potentials. Based on this observation, it is argued that the Local Variational Principle 
performs better than the equivalent methods based on the centroid path idea and on the Gibbs- 
Bogoliubov-Feynman variational principle, especially in the range of low temperatures. 
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I. INTRODUCTION 

The Gibbs-Bogoliubov-Feynman inequality (GBF) is a 
restatement of the second law of thermodynamics. How- 
ever, the motivation of the present work is the equally 
important fact that the inequality provides a variational 
approximation to the Helmholtz free energy. Historically, 
Gibbs first stated the inequality for classical systems, 
then Bogoliubov and Feynman generalized it to quantum 
systems in the operator and the path-integral formalism 
of quantum mechanics, respectively. Perhaps at the ex- 
pense of losing the original physical significance, the local 
variational principle I develop in this work is intended to 
be a mathematical basis for the design of more efficient 
computational methods dealing with statistical quantum 
systems, with special concern for their low temperature 
behavior. To fully justify the need for a local principle, 
we first have to give a short review of the Feynman and 
Gibbs-Bogoliubov inequalities and at this point, we shall 
also introduce some notations of use throughout the pa- 
per. 

The path integral formulation of the statistical me- 
chanics began with the Feynman's realization at an "intu- 
itive" level that the density matrix of a monodimcnsional 
quantum particle is the expectation value of a suitable 
function of a Brownian Motion |lj . Feynman was actually 
working on the real time Schrodinger equation, but for 
the imaginary time analog the theory was made rigorous 
by Kag |2|, the product being the well known Feynman- 
Kag representation formula (Theorem 6.6 of Ref. |3J) 



p{x,x';(3) 
Pf ree (x,x';f3) 

where 



,exp ■ 
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x(t) 



dt 



(1) 



and 



p free (x,x';0) - y/ 



exp 



(x-x'f 



is the density matrix for a similar free particle. B® de- 
notes a standard Brownian Bridge (see p. 40-41 of Ref. 3 
and p. 430-431 of Ref. 0) and the expected value in (JJJ is 
taken with respect to its underlying probability measure. 

For the sake of simplicity, we shall be concerned mainly 
with the monodimcnsional case, but the reader should 
observe that the theory is in no way restricted to this 
case. This is so because the Feynman-Kag formula has 
a straightforward multidimensional generalization: one 
simply utilizes an independent Brownian Bridge for each 
physical degree of freedom. However, we explicitly ad- 
dress various multidimensional problems, whenever they 
significantly differ from their monodimensional version. 
As stated, the main theorems obtained in this paper re- 
main true for the multidimensional systems. 

The Fourier Path Integral (FPI) implementation of Q , 
which we exclusively use in this paper, is due to Doll and 
Freeman |5j and is based on the exact representation of 
the Brownian Bridge as a Random Fourier series with 
the coefficients being independent identically distributed 
(i.i.d.) Gaussian variables. To rephrase their result in 
the spirit of the Feynman-Kag representation formula, if 
fl is the space of infinite sequences a = {a\, a.2, ■ ■ ■ ) and 



P[a] = J] /iK) 

fc=i 



(2) 



is the (unique) probability measure on fl such that the 
coordinate maps a — » a k are i.i.d. variables with distri- 
bution probability 



H(a k € A) 



1 



v(t) =x +(x' - x)t 



/2tt 



- 2 / 2dz 



(3) 



then, 
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is equal in distribution to a standard Brownian Bridge. 
Let us introduce the path-averaged potential functional 



U(x, x 



where 



„1 oo 

',a;/3)= / V[x(t) + J2a k a k sm(kiTt)]dt, (5) 
Jo fe=l 



2/3fi 2 1 



7r 2 m fc 2 ' 

and make the convention that whenever x = x', the prime 
a; is dropped so that U [x, a; 0) = U (x, x, a; (3). With this 
notation, the FPI version of the Feynman-Kag represen- 
tation formula is 

p{x,x';(3) = p free {x,x';f3) 

x / dP[a]exp[-(3U(x,x',a;(3)} . (6) 

In his treatment of the Frohlich polaron problem 
Feynman constructed an upper bound to the free energy 
of a quantum system by means of the inequality (see 
formulae 3.52 and 3.53 in Ref. 0) 



F<F{ + {U-U{) ULh 
where, in general, (0)s>_ stands for the average 



(0)i 



J R dxJ n dP[a]e-^^0(x,d;l3) 



f K dx J n dP{a]e 



-0UUx,a;l3) 



(7) 



(8) 



The functional Ur(x,x' , a; /3) was taken to be of the 
form (J5J) for some trial potential VJ(x) depending upon 

a set of parameters b = (b±, 62, ■ ■ ■)> but this is not a re- 
quirement and essentially any function satisfying some 
mild integrability conditions can be utilized in the Feyn- 
man inequality Q. 

As argued by Feynman (see Chapter 11 in Ref. 8), the 
zero temperature limit of the inequality Q is 



e < 



<P° b \H\<P° b ) 

mm ' 



(9) 



where 0° is the groundstate eigenfunction of the trial 
Hamiltonian 



h 2 d 2 



(10) 



Eq. [5] is in agreement with the Raylcigh-Ritz principle 
for groundstate eigenfunctions and, as an approxima- 
tion, can be arbitrarily sharpened by use of more accu- 
rate trial potentials. By the inherent continuity of such 
problems, these good variational estimates of the ground- 
state energy imply good estimates of the Helmholtz free 
energy for the entire low temperature regime, fact hard 
to achieve by other means. This helps explain the suc- 
cessful application of the Feynman variational princi- 
ple in a variety of theories dealing with the evalua- 
tion of the thermodynamic properties of quantum sys- 
tems 0, E3, EH E3, EE E3 



The Gibbs-Bogoliubov inequality [1 
following bound to the free energy 



16] provides the 



F<F-l 



Tr\{H - Hfie-rfi] 

a £11 ' 



Tr(e~ pA i) 



(11) 



which for spinless particles is proven to be equal to the 
one given by the Feynman inequality Q, whenever the 
functional U((x,x' ,a; j3) can be cast in the form of the 
equation (J5J for a given trial potential Vr(x). In this 
situation one talks about the Gibbs-Bogoliubov-Feynman 
inequality ( GBF) and of the corresponding variational 
principle consisting of the minimization of the right- 
hand expressions in the formulae Q and (|ll|l on the 
set of parameters b. The reader should not conclude that 
the Gibbs-Bogoliubov inequality is automatically weaker 
than its path-integral counterpart. For instance, in the 
case of a fermionic system, the Gibbs-Bogoliubov in- 
equality is still true if the trace is restricted to the Hilbcrt 
space of antisymmetric functions, with the slight require- 
ment that the trial Hamiltonian be totally symmet- 
rical under the permutation of identical particles. How- 
ever, there is no known path-integral equivalent to the 
resulting inequality, the difficulty being related to the 
so-called fermionic sign problem [l7|. It is for this rea- 
son that we shall restrict the development of our local 
variational principle to spinless systems. 

The GBF usefulness depends upon our ability to ana- 
lytically compute the integrals on the right-hand side of 
the equation J7J, at least the ones with respect to the 
Fourier coefficients. This effectively restricts the choice 
of trial potentials to a handful (in most cases a quadratic 
potential) and it is in poor match with the fact that the 
Feynman estimate is global, involving an integration over 
the physical coordinates. We thus arrive at the first mo- 
tivation for our work: a local fitting (pointwise in the 
configuration space), as opposed to a global one, would 
make more use of a simple reference potential. Then, a 
pointwise approximation of the density matrix can always 
be improved by other means than the use of more compli- 
cated trial potentials, the default choice being the Trotter 
composition rule Finally, we will show that the lo- 
cal variational principle developed in this work provides 
new information about the density matrix, information 
unattainable from GBF. Our perspective on the compu- 
tation of the density matrix is thus changed: instead of 
seeking better reference potentials, we try to find the 
variational principle which makes the best use of a given 
reference potential. 



II. THE LOCAL VARIATIONAL PRINCIPLE AS 
A LOWER BOUND FOR THE FINITE 
TEMPERATURE DENSITY MATRIX 

The purpose of this section is to define the Local Varia- 
tional Principle (LVP) and further justify its importance. 
In addition, we shall consider the particular case of LVP 



3 



when the reference potential is the quadratic one and 
compare this case with the centroid based approxima- 
tions 

[3 H ED, Elm El , p articularl y with EFLT & 

which is the GBF analog. 

The Gibbs-Bogoliubov-Feynman inequality is a conse- 
quence of Jensen's inequality and I remind the reader the 
latter's statement (see p. 14 in Ref. 0): 

Theorem 1 (Jensen's inequality) If(fl,P) is a prob- 
ability space, if g : Q — > (a, b) is integrable and if F is 
convex on (a,b) with — oo < a < b < oo, then 



Fog dP > F gdP 



By default, whenever wc apply Jensen's inequality in 
this work, it is understood that the convex function is 
the exponential F(x) = exp(— x). 

To begin with the definition of the Local Variational 
Principle, let us perform a change of measure in JSJ of 
the form 

p(x,x';P)=p free (x,x';(3) [ dP[d} e -^ u ^ x '^ 



x exp{-/3 [U(x, x', a; (3) - U&x, x' , a; p)] }, (12) 

where UUx, x', a; 0) is any measurable function depend- 
ing upon a set of parameters b — (&i, 6 2 , . . .) such that 



p' l (x,x';(3)=p free (x,x';(3) / dP[a} e -^ u ^ x '^ (13) 



has an integrable diagonal. Defining a new probability 
measure by the relation 



dP! , r a\ fa 



Pfreejx, X'] j3) c _g^( a; ,x , ,a;ffl d p[ 5 ] i ^ 



p'Ax,x';P) 



we may rewrite (|12ll as 



p{x,x';(3) = pi{x,x'\(3) / dP( xx ,- b 



x exp{-/3 [U{x,x',d;/3) - U L b (x, x' , a; j3)\ } . (15) 

Now, use of the Jensen's inequality produces the local 
variational inequality 



p{x,x';0)>pl{x,x';f3), 



(16) 



where 



pf{x,x';[3) = pi(x,x';0)exp{ - f3 J dP ( ' x ^^[a] 

x[U(x,x',a;p)-U-l(x,x',a;p)]}. (17) 

This inequality is the cornerstone of the variational meth- 
ods, providing a bound from below to the density ma- 
trix. Specialized versions of the inequality were consid- 
ered before whether as a starting point for the definition 



of the Partial Averaging method [24( or in the context of 
the Fey nman-Kleinert variational-perturbational theory 
[25L l2fij . We should remark here that the nonnegativ- 
ity of the density matrix, which stems from the reality 
of the path-averaged potential functional, played an im- 
portant role. Therefore, the inequality is not true for 
general complex /3. The local variational inequality i|ltj|) 
implies the Feynman inequality (JJJ. The latter can be 
deduced by setting x = x' in 1|16|) . integrating over x, 
working along the same lines as in (|14I15I) . and finally 
using again Jensen's inequality to obtain 



e -0F > e - P Fi 



{-" 

[U(x,a;p)-Ui(x,a;p)} } 



> e 



-0H 



exp 



(18) 



which produces Q upon taking the logarithm. 

For the rest of the paper, we shall only be concerned 
with the case when the functional Ul(x, x' , a; (3) is the 
path average of some reference potential depending upon 
the set of parameters b 



Ur(x,x' \a;(3) 



„i 00 

/ V{[x(t)+Y^a k a k sm(kirt)]dt. (19) 
Jo k=i 



The reference potential is assumed to be a bound- 
ing potential, with a discrete spectrum and a unique 
and strictly positive groundstate eigenfunction. We shall 
denote its eigenfunctions by <t>^(z) and the correspond- 
ing eigenvalues by e^'. Taking the supremum in the 
equation (|16fl over the set of parameters b produces the 
sharper Local Variational Principle (LVP) 

p(x,x';f3) > p best (x,x';[3) = sup p^(x, x'\ 13). (20) 

6 

We say that pb es t(x,x'; (3) is the best approximation of 
the density matrix in the sense of LVP. The extremum of 
the maximization problem l|20() is attained on some pa- 
rameters b = B{x,x' , (3) which generally are functions of 
position and temperature (in case of multiple maxima, 
choose arbitrarily one of them) and as a direct conse- 
quence, pbest(x , x' '; (3) is no longer the density matrix of 
a trial potential. 

At this point, it is useful to consider the special case 
of LVP when the reference potential is the harmonic os- 
cillator one. The method will be termed HO-LVP. I only 
present the monodimensional version as a clear sugges- 
tion of how the inequality H2(J|I can be employed. The 
multidimensional version as well as the specific numeri- 
cal implementation and the related problems will make 
the object of a separate paper. The most general monodi- 
mensional quadratic potential has the form: 



V^x) = -m u; 2 (x-z) 2 



(21) 
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where the translational variable z and the frequency u> by means of the FPI path integral formulation, give the 
are the parameters to be determined by the local varia- following HO-LVP approximation for the density matrix: 
tional principle. Straightforward but lengthy calculations 



Pt ^ X ''f = ho (pC) exp f - IfPlPhBtfC) - lp 3 A 2 h 6 ((3C) + l/3 2 C%^C) 

Pfree\X, X iP) L Z Z Z 



x exp ■ 



(3 / df Vt !U1 x{t) 



2h 2 



f3 2 Bh 2 {(3C, t) 



2fr? 



(3 2 Ah 3 (f3C,t) 



(22) 



r 



In the above, 



C = - 

IT 



B = 
A = 



2W- 



-(x + x' - 2z) 
'-(x-x'). 



(23) 



Even though A, B and C are in fact functions of u>, z, 
x and x', we do not write their arguments explicitly in 
order to save typographical space. The h-functions are 
tabled in Appendix A. In addition, 



VtA x ) = / V(x + z) 



1 



exp 



2T£(t) 



dz 



(24) 

is a convolution of the original potential with a Gaussian 
of width 



r£(t) 



(25) 



and is called a Gaussian transform of the potential. 
The inequality 12U|) simply states that 

p(x,x , ;(3)>p a z Jx,x';(3) 

so that, for each pair of points (x, x') the maximum 
of the right-hand side expression is attained on some 
optimum values of the parameters z = Z(x,x';/3) and 
w = Q{x,x';j3). 

I do not discuss here how the HO-LVP technique can 
be numerically implemented for practical applications. 
With the sole difference that there are monodimensional 
integrals against the parameter t to be computed numer- 
ically (a tractable problem), the HO-LVP is on par as 
computational difficulty with the centroid based meth- 
ods [2£j, l2ll I22L l23l ] and it is amenable to similar approx- 
imations (see Ref. ^3 for an example) . They all involve 
local minimizations and integrations in the configuration 
space, respectively in the centroid space. 

Rather, we shall emphasize the differences between 
such methods. Since it is required that the density ma- 
trix of the reference potential be analytically known, a 



common feature of the methods is the fact that only 
simple references, as for instance the quadratic poten- 
tial reference, are computationally feasible. Therefore, 
it is desirable that the approximation which makes the 
better use of the simple reference potential be employed 
in actual simulations. As such, HO-LVP has two impor- 
tant advantages over EFLT [2J: (a) it can be arbitrarily 
improved by Trotter composition |18j | and (b) it gives an 
approximation to the true density matrix which provides 
more information about the system than the variational 
approximation to the centroid density matrix. 

Because the first property is clear, we shall be mainly 
concerned in this paper with proving the second asser- 
tion. To this point, we notice that the high temperature 
limit of any density matrix is the classical one. However, 
as the temperature is lowered, the discrepancy between 
the classical and the quantum density matrices increases 
because the thermodynamic spread of the paths entering 
the Feynman-Kac formula also increases. In fact, if the 
number of variables used to parameterize the paths is 
kept constant, the thermodynamic energy estimator for 
non- variational methods as DPI [H ^ El H13 or FPI 
extrapolates to the classical energy in the low tempera- 
ture limit too [H]] i a phenomenon dubbed "classical col- 
lapse." Giachetti and Tognetti as well as Feynman 
and Kleinert ^(J noticed that this is not true of the varia- 
tional methods based upon the GBF principle. Following 
their line of thought, one may argue that the low temper- 
ature limit of the energy estimator for the EFLT centroid 
method (see the equations 2.34, 2.35, 2.41, and 2.42 of 
Ref. ED) is 



inf ■ 



CI 



> eo, 



i.e, the expected energy of the best Gaussian wavepacket. 
[Remember that <jfi z w (x) is the groundstate eigenfunction 
of the reference potential given by Eq. |^. This is a 
quite remarkable fact because the harmonic oscillator is 
known to be a good approximation of the potential sur- 
face around the main local minima. 

However, one very important aspect of the centroid 
density matrix p c (x;/3) is that it bears no direct con- 
nection to the true density matrix. For instance, given 
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p c (x; j3), one cannot compute exactly the ensemble aver- 
age potential energy 

J R p(x-/3)dx 

though useful approximations are known |22j |. As far as 
the total energy is concerned, this can be exactly evalu- 
ated with the help of the T-method estimator (see Sec- 
tion V for definition). But up to some functionals of it as 
for example the partition function, this is the only prop- 
erty that can be computed exactly once the centroid den- 
sity matrix is known. Clearly, there is no H-method esti- 
mator for the centroid density matrix. Things are totally 
different in the case of the Local Variational Principle be- 
cause this provides an approximation for the true density 
matrix and so, the expectation values of different opera- 
tors are readily available. Even more, we shall later show 
that the zero temperature limit of the H-method estima- 
tor is a groundstate energy estimate always better than 
the corresponding centroid one, the latter being matched 
by the low temperature limit of the LVP T-estimator. 
I hope this would be enough evidence to convince the 
reader that LVP provides the better description of the 
physical system. 

From the above discussion, we infer that the quality of 
a variational approximation is dictated by its low tem- 
perature limit, and in the next section we shall estab- 
lish what this limit is in the case of the Local Varia- 
tional Principle. The reader should not forget that LVP 
provides a variational bound from below to the finite- 
temperature density matrix and that it is intended as an 
approximation method for this density matrix. There- 
fore, LVP is in no way limited to the computation of the 
groundstate eigenfunction, which is however our object 
of interest for the next section. 



III. THE EIGENFUNCTION 
REPRESENTATION OF THE LVP 

In order to establish the low temperature limit of the 
Local Variational Principle, as well as to show that the 
local variational inequality (|2U[1 can also be interpreted as 
a generalization of the Gibbs-Bogoliubov inequality, we 
need to express p^(x, x'; /3) in terms of the eigenfunctions 
and the eigenvalues of the potential V{ (x). We again de- 
velop the theory in full generality, rather than discussing 
the special HO-LVP case. 

In anticipation of the final result, let us see that 
the eigenfunctions and the eigenvalues of the perturbed 
Hamiltonian 



respectively, 



H{ + \ [V(x) - V{(x)] 



(26) 



as given by the first order Rayleigh-Schrodinger pertur- 
bation theory are of the form 



(27) 



' ~ f * + X($\V-V{\$), 



(28) 



provided that the eigenfunctions <t>\{x) are chosen such 
that the perturbation V(x) — Vt(x) is diagonalized on 
each degenerate subspace. Letting D k = {i G N : = 
e^} and = {4>i\V — V^'|<^), the exact expressions for 



the coefficients Cki are (see Chapter 5 in Ref. 32): 



Cki 



^ki 



e i _ gfc 

b b 



4*4 



Vn - Vkk ^ d - s? 



v kk ~~Z Ej — el 

3 <£D k b b 



Cki = 



, e\ = £,v u = Vkk- (29) 



I warn the reader that the equations l|27l) and (|28|l are 
exact to the first order in A, for instance, 



e b,\ e b I .fci T / \rl\ ±k\ 



lim 

a^o A 



With these preparations, we are ready to prove an im- 
portant lemma: 

Lemma 1 



P J dP ( x ,*',b,P) ^ U ( x ' x', a; (3) - Ufa, x' , a; /?)] 



e-^l 



EZo 4(x)^(xQ(^|V - V{\^ k b )e-H 



(30) 



Proof: Write V{(x) = V(x) - VI {x) and 

Ut'(x, x' , a; (3) = U(x, x' , a; ff) — Ul(x, x' , a; f3), 
so that: 

„1 oo 

U%(x, x',a;{3) = / V{ [x(t) + V a k a k sin(fcirf)] dt. 
Jo fe=i 

Then consider the equality: 

1-exp [-\/3U%{x,x',a;l3)] 



(3U%{x,x',a;(3) = lim 



A 



Remembering the definition l|14|) of the probability mea- 
sure Pi x , i.0^fl\ and the eigenfunction series representa- 
tion of a density matrix, we learn that 
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where <fi^ x {z) and A are the eigenfunctions, respectively 

the eigenvalues of the perturbed Hamiltonian l|26|l . 

For small A, we may use the Rayleigh-Schrodinger per- 
turbation theory to compute the spectrum of the per- 
turbed Hamiltonian H{ . . The reader should realize that 

0,A 

the only corrections needed are the ones to the first order 
in A (which are exactly given by the Rayleigh-Schrodinger 
perturbation theory), since the others will cancel upon 
letting A go to zero. To conclude the proof, use the for- 
mulae <|27I28[1 and explicitly compute the limit in 131|) . □ 
Lemma 1 together with the well-known series repre- 
sentation of a density matrix 



p> h (x,x>;f3) =J2$(x)$(x>)e^ 



(32) 



k=0 



essentially solves the eigenfunction representation prob- 
lem. The functions 



and 



-0e- 

e p <• 



(33) 



if {x 1 x' ] p)=Y J ^(x)^{x'){cj ) l\V -V{\4>l)e 



04 



k=0 



have the following obvious properties: 



T^\x;f3)dx = 



(34) 



(35) 



and respectively 

p CO 

/ T^(x;(3)dx = ^l\ V -VM)e-^ 
Jr z n 



fe=0 



(36) 



Therefore, 

Theorem 2 The eigenfunction expansion form of the lo- 
cal variational inequality \2U\) is 

p(x,x';(3) > pf(x,x , ;l3)=p! B (x,x';0) 
x exp 



T&(x,x';l3)+prP(x,x?;0) 



•(37) 



Before continuing, the reader is advised to ponder over 
the value of this theorem by analyzing the HO-LVP ap- 
proximation p" u (x, x'\ (3) compactly given as a monodi- 
mensional integral against t by Eq. 1221 The same HO- 
LVP approximation can be exactly written in terms of 



the eigenvalues 



and the eigenfunctions 



^(») = (2 fc *!)- 1 /»(^) 1/ V/»^(0 

of the harmonic oscillator reference potential l|21|l , in the 
form given by Theorem 2. Here, £ = [m^w jK) x ^ 2 (x — z), 
while Hk(x) stands for the respective Hermite polyno- 
mial. Of course, the eigenfunction representation is of no 
practical use, but it allows us to study the low tempera- 
ture behavior of the density matrix. 

In the form l|37|l . the local variational inequality is 
readily seen to imply the Gibbs-Bogoliubov inequality. 
Indeed, setting x — x', integrating over x and applying 
Jensen's inequality produces: 



-PF > e -PH / dx 



-0F-L 



x exp 



T^(x;P)+(3Tf\x;l3) 



> e~ pF i exp I -13- 



Tr 



Tr(e- pH i) 



, (38) 



where we used the relations (|35|) and l|36(l . The last equa- 
tion produces the Gibbs-Bogoliubov inequality upon tak- 
ing the logarithm. Moreover, since we performed the 
same operations as for <|18ll . we also get a proof of the 
equivalence between Feynman and Gibbs-Bogoliubov in- 
equalities, provided that the form (|19|l for the trial poten- 
tial is assumed. It is in this respect that we regard LVP 
as a generalization of both aforementioned inequalities, 
even if the best density matrix predicted is not necessar- 
ily derivable from a potential. 

The remainder of this section deals with the low 
temperature behavior of the LVP density matrix 
Pbest{x,x'; (3). An immediate corollary of Lemma 1 is 
the equality 



lim 



U (x, x' , a; 13) 



-Ui(x,x',a;{3)\ -(4\V-V{\4) 
-- S- b (x) + S- b (x% CW) 
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where 



fe=i 



4 



4 



(40) 



is a function which does not depend upon temperature. 
In deducing l|40|) . one uses the fact that the groundstate 
eigenfunction of the trial potential V b (x) is not degener- 
ate. Then, the asymptotic formula 

implies 

pl(x,x';P) « 4(x)4>l(x')exp{-[S- b (x)+S- b (x')]} 

xexp[-/3(eg + (^|F-^|^))].(41) 

I warn the reader that here and in the remainder of the 
paper the sign w is used to denote a low temperature 
asymptotic form, its rigorous interpretation being: 

Urn {exp [(3 (e° b + ($\V - V{\$))] p%(x,x'; (3)} = 

= <t>l(x)4(x')e X p{-[S- b (x)+S- b (x')]}. 

Looking at the equation 1)41(1 . we see that the factor 
containing the S b (x) functions simply disappears in the 
original GBF equation because of the identity ()35|1 . Thus, 
our theory brings some new information about the shape 
of the groundstate density matrix, and we shall later 
prove that the correction factor is always an improve- 
ment in the energetic sense. After an obvious simplifica- 
tion of the terms explicitly involving the potential V4 (x), 
the following theorem is immediate: 

Theorem 3 The asymptotic formula of p best {x, x'; j3) at 
low temperature is: 

p best (x,x';(3) « ^(x)^(x')exp{-[5 6 (x) + ^(^)]} 



x exp [-(iE{4>i)] 



b—B(x,x',oo) 



where 
E(tl>) 



— \\V^x)\\ 2 +^(x) 2 V(x) 



dx 



(42) 



(43) 



and the functions B{x, x', oo) are computed by the follow- 
ing recipe: 

1. Minimize the functional E (ip^) . If it is unique, the 
value of b on which the minimum is attained 
becomes B(x, x , oo) Vx, x' € M. 

2. If there are multiple minima of E{4>^), pick an 
arbitrary one among those that further maximizes 

4(x)4(x')e X p{-[S b (x)+S b (x')]} 

at each pair of points (x,x'). 



Let us analyze a little more closely what the last the- 
orem says. Assume we are in the simple case when the 
minimum of the functional -E(0j) is unique. Up to a 
normalization factor, Theorem 3 predicts the following 
approximation to the groundstate eigenfunction: 



^b{x) = 4>i (x) exp 



$(x) (<f£\V-V{m 



k=l > 



4 



(44) 

where the optimal parameters b do not depend upon the 
coordinates (x, x'). Thus, the Theorem 3 does not simply 
predict the function 4>^{x), though the thermodynamic 
weight is computed with respect to this function. An im- 
mediate question is in place: What can we say about the 
quality of the above eigenfunction? The quite remarkable 
answer is proved in the next section (see Eq. 1 79(1 , and 
says that the expected energy of ip b (x) is always smaller 
or equal to the expected energy of 4> b {x). In other words, 
LVP predicts an energetically better groundstate eigen- 
function, and we shall prove in Section V that we can 
recover its expected energy by use of the H-estimator. Fi- 
nally, for multidimensional systems, LVP predicts a cor- 
related approximation of the groundstate eigenfunction 
even if the reference is a sum of monoparticle potentials. 
The low temperature density matrix given by Ea. 1421 has 
no GBF equivalent and justifies our claim that LVP is a 
separate and more powerful principle. 



IV. THE EXPECTED ENERGY OF THE LVP 
GROUNDSTATE DENSITY MATRIX 

Let us remember that our special interest for the 
groundstate density matrix is due to our experience that 
the various approximations used to compute the finite- 
temperature statistical properties of a physical system 
worsen in the low temperature regime. By the intrin- 
sic continuity of the variational methods (see Section V 
for further clarifications), a good approximation of the 
groundstate density matrix necessarily implies a good ap- 
proximation for the finite-temperature density matrix. In 
this section, we shall analyze the expected energy of the 
groundstate density matrix predicted by Theorem 3, but 
we assume a special form of the density matrix which is 
encountered in practical applications whenever the po- 
tential V{x) has a finite number of local minima. 

There is one special parameter bo which accounts for 
a translation and which we add to the list of parame- 
ters b = (61,62, .. .)• From now on, we shall conform to 
the convention that if not written explicitly in an expres- 
sion, 60 is assumed to be part of the list of parameters 
6 i.e., 6 = (60, 61 . . .). Otherwise, if 60 does appear in an 
expression, the list 6 is assumed not to contain it. The 
importance of this parameter consists of the fact that, 
if it is included, the optimizing coefficients B(x, x'; 00) 
usually become constant on certain regions of the config- 
uration space, which are identified with the main wells 
of the potential. Of course, for an n-dimensional sys- 
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tern, there are n translational parameters, one for each 
dimension. 

To begin with, we replace (|19|l by the slightly more 
general form 



U' l (x,x',a;/3) = / V^[-b Q +x{t)+J2 &k°k sin(fc7rt)] di. 
Jo fc=i 

(45) 

According to our convention, V~(x) does not depend ex- 
plicitly upon bo, the value of this parameter, which sets 
the origin of the potential, being automatically deter- 
mined by LVP. All our results remain true if computed 
with respect to the local reference potential V£ b (x) = 

V£(x — bo). If the eigenfunctions of ^(x) are </>^ ^(x) 



and the corresponding eigenvalues are ^ then we have 
to introduce the two local quantities: 



V z (x) = V(x + z) (46) 



and 



^1 
2m 1 



Vip(x)\\ 2 + ip(x) 2 V(x + z) dx. (47) 



Also, we shall use z instead of bo and let A index all the 
pairs (z a , B a ) on which the minimum of the problem 



E best =mfE z (4) 

z,b 



(48) 



is achieved. The new system of indexation makes the 
old one superfluous, so we shall drop some indices. We 
define: 



S a (x) 



<t> k a (x-z a ) (4>l\V, a -V' B \<fc) 



^ 4>° a (x - Z a ) 



and 



ip a (x) = 0° (x - z a ) exp[-5 Q (x)], 



(49) 



(50) 



where (fraix) and e a are the eigenfunctions, respectively 
the eigenvalues of the trial potential Vg (x). 

With these notations, Theorem 3 takes on the special 
form: 

Theorem 4 // \4-5\j is assumed, then the asymptotic for- 
mula of pbest(x,x'; (3) at low temperature is: 

p best (x,x';/3) w exp(-(3E best )pl est (x,x'), (51) 

where E best is the defined by and 

Pbest( x ,x') = SU P 1pa(x)^ a (x'). 



(52) 



To appreciate the importance of the translational pa- 
rameter z, let us perform the minimization in (|48|l in two 
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x(au) 

FIG. 1: A plot of the symmetric double well quartic poten- 
tial 1561 and of its associated effective potential 1531 . 



separate steps. Firstly, we construct an effective poten- 
tial 



V ef (z)=ME z (4) 

b 

and secondly, we compute 

E best = inf Vef(z). 



(53) 



(54) 



For monodimensional systems, it is usually the case that 
the minimum of the first problem is attained on unique 
points b = B{z) while for multidimensional ones (espe- 
cially for systems in condensed phase) there is usually 
a finite number of minimizing parameters. The effective 
potential is in fact a mollification of the original poten- 
tial V(z), to which it converges as the ratio h 2 /m goes to 
zero. For systems in condensed phase, it is often the case 
that both the original and the effective potentials have 
finitely many global minima, and in this section we shall 
assume that there are finitely many pairs (zi,Bi),i £ 1, TV 
on which mf z b E z (<fi®) is attained. A more general result 
will be proved in Section V. If we set 

A = {{x,x') e R 2 : p° best (x,x') = Mx)M*')} , (55) 

then the sets Di are assumed to be disjoint except for 
their (topological) frontiers which are required to have 
measure zero. It follows that the optimizing coefficients 
are constant on the interior of the sets Di and that 
p best (x,x') is twice derivable with continuous derivatives 
on the same interiors, yet continuous on the entire plane 
M 2 . Therefore, the diagonal density matrix p best (x) as 
well as its square root have similar continuity properties 
with respect to the diagonal sets Df, defined as the in- 
tersections of the Di's with the line of equation x = x' 
[remember the convention p best {x) = p best (x, x)]. 

The rest of this section deals with the evaluation of the 
expected energy of p best {x, x'). To reinforce the proofs, 
we study a simple example of a quartic double well poten- 
tial in the context of the HO-LVP approximation, along 
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with the general approach. We shall set h — 1 , and con- 
sider a particle of mass m = 1 moving in the potential 



V(x) = -{x- A) 2 {x + A) 2 , 



(56) 



where A — 1.5 (see Fig. 1). The reference potential is a 
quadratic one, of variable frequency to > 0: 



1 



V^(x) = ^mu 2 x 2 . (57) 



The functional (|47l) can be worked out explicitly to be 

(58) 



_ . , u Zz 2 -A 2 



_3 



Fig. 1 also contains a plot of the effective potential 
V ef (z)= ME t (u) 

ui>0 

and shows that V e f{z) attains its global minimum 
Ei> es t = 1.404 on the two symmetric points Z\ = —1.292 
and Z2 — 1.292. The corresponding optimum reference 
potential frequencies are lo\ — ll>2 = u> = 2.584, equal by 
the symmetry of the problem. The Si(x) functions can 
be expressed in terms of the Hermite polynomials as 



fc=3,4 



-0.539, c 2 = 0.539 and c\ = <% = 0.092. In 



where c\ — 

general, it can be shown that all the coefficients Ck with 
k > 2n vanish for any polynomial potential of rank at 
most 2n, while the coefficients c\ and C2 vanish for all 
potentials. The functions ipi(x) have the form 

ipi(x) — exp -~w(x - Zi) 2 - Si(x) , (60) 

so that 

Pbest( x ^')= max ^,(1)^(1'). (61) 
*e{i.2} 

In Fig. 2, one may see that the low temperature density 
matrix predicted by LVP is symmetrical at reflection with 
respect to both the main and the secondary axis. Though 
continuous on the entire plane, the density matrix has a 
cusp along the secondary axis. The sets Di are readily 
identified: D\ = {(x,x') : x' < x} and D 2 = {(x, x') : 
x 1 > x} with the diagonals D\ = {x < 0} and -DJ = 
{x > 0}. 

Let us go back to the energy evaluation problem. The 
lack of continuity of the first derivatives on the bound- 
aries dDi requires a careful analysis of the kinetic energy. 
We consider two estimators: 



Ki{p) = 



h 2 k 



dxdx 



-P(x,x')\ dx 



2m 



J R p(x)dx 



(62) 



\ 
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FIG. 2: A plot of the low temperature density matrix pre- 
dicted by LVP. There are only two maxima instead of four 
symmetrical ones, the true density matrix would present. 



and 



K 2 (p) 



h 2 L 7r^zp(x, x')\ ,dx 

11 JR dx A r\ 1 I I x =x' 

2m J R p{x)dx 



(63) 



where the derivatives are regarded as functions almost 
everywhere defined and not as distributions. We shall 
denote by E\ (p) and E 2 (p) the associated energy estima- 
tors, obtained by adding the expected potential energy. 

A little thought shows that the following equality holds 
for all the points (x,x') in the first and the third quad- 
rants of the plane: 

pLA^ = y/plstWy/plstW- ( 64 ) 

In fact, the square root of the diagonal density p1 est {x), 
which has a cusp in the origin as shown in Fig. 3, will 
play such an important role in our development that it 




FIG. 3: For the double well quartic potential (15611 . the ap- 
proximate groundstate wavefunction tp (x) defined by 16511 
has a cusp at the origin. 



10 



deserves a notation: 



^o(x) = ^P° best (x). 



(65) 



In general, one may show that around each diagonal 
point (x,x) on the interior of some set Df, there is a 
small neighborhood in R 2 , say the ball B[e, (x,x)], such 
that (|64|l holds on B[e,(x,x)]. Relation 1)64(1 needs not 
hold for the diagonal points (x, x) that are precisely on 
some frontier dDf but from the point of view of integra- 
tion theory, this does not matter because the frontier has 
measure zero. Consequently, the following equalities are 
true: 



and 



Ki(f° best ) 



K l{Pbest) 



H 2 J R ||V-0 o (a;)|| dx 
2m J R ip (x) 2 dx 



h 2 J R ip {x)-£jip (x)dx 
2m f R ip (x) 2 dx 



(66) 



(67) 



Because p best (x) is continuous on the entire R by the 
way it was constructed, it can be proven that ip (x) is in 
the Sobolev space iJ 1,2 (R), thus a permissible trial func- 
tion for the groundstate eigcnfunction of the potential 
V(x). For functions ip{x) in iJ x ' 2 (R), the Rayleigh-Ritz 
principle states that 

h 2 J R \\^(x)\\ 2 dx J R V(x)iP(x) 2 dx 

e < r i , no i 1 r i , no i ■ (0») 



2m f R ip(x) 2 dx 



f R ip(x) 2 dx 



Consequently, the correct variational definition of the 
kinetic energy is given by the formula l|62[l and we have 
our first important result: 



e < E x (p° beat ). 



(69) 



For the case of the quartic potential, the exact values 
are: e = 1.292 and Ex{p best ) = 1.342. We see that our 
estimation of the groundstate energy is better than the 
one given by GBF, which is E best = 1.404. We shall 
prove that this is no mistake, and that the energy of 
the asymptotical low temperature density matrix pre- 
dicted by LVP is always lower than the one predicted 
by the GBF. We do this in two steps: first we prove that 
Ei(p° best ) < E 2 (p° best ) and then that E 2 (p° best ) < E best . 
Integration by parts produces 



d 2 

^o{x)-^tlJo(x)dx 



2 ^ 



dx 



|VV>o(x)|| dx, 



(Xi) 



(70) 



where the points Xi are separating two consecutive sets 
Df and Df +1 on which the diagonal density matrix takes 
the values ipi( x ) 2 an d tpi+i(x) 2 respectively. Notice that 
ipi{xi) 2 = il) i+ i(xi) 2 and that 

if>i+i(xi + h) 2 > ipi(xi + h) 2 



for all positive and small enough h or, otherwise, 
p best (x) = ipi(x) 2 for some x G Df +1 contradicting the 
definition of the set. Therefore, 



dx 



(Xi) 



lim 

h\0 



i/>? +1 (xi + h)-Tp?(xi + h) 
h 



> 0, 
(71) 

(72) 



which together with (|70fl proves that 

Ei{p° b est) < E 2 {p best ). 

For multidimensional systems, the same reasoning can be 
performed along the normals Pjj to the surfaces dDf n 
dDJ separating the sets Df and DJ. The normal i/y is 
assumed oriented from the Df toward the Df set. Then, 
the analog of (j7I|) is 

(73) 



V [V>] - V?] • Vii > on dDf n dDJ 



and the analog of (|70|l is 

ipo(x)Aip (x)dx — / || V-0 o (a;)|| 2 da; 

proving again J75J. 

Finally, let us show that E 2 (p best ) < E best . Because <^>° 
is strictly positive, we can write any other eigcnfunction 
4>i(x) as the product ff{x)4>®(x). By direct substitution, 
one can show that the function f^(x) satisfies the follow- 
ing equation: 

4 A/ ' (l) 4 V ^ (l)Jl ' V/ ' W = (ti-z°i)fi(x)- 

(75) 

It follows then that 



2m 



ASi(x) - —Vh#°(x- z t ) 2 } ■ VSi(x) 



2m 



E 



f^ x - Zi )(^\V Zi -Vg\^). (76) 



The sum of the last series equals 

V(x)-V^(x-z i )-(cl> i \V z -V^) 

by the completeness of the system of eigenfunctions 
{<fii(x);k > 0} and the translational invariance of the 
integrals involved, so we end up with the equality 



^(x)~Vh[^ 
= [V(x)-VUx- Zi )]-(4 



-z^-VSiix) 



(77) 



With the help of (|77|l . one can show by explicit compu- 
tation that 

h 2 

-Aik(x) +V{x)ipi{x) 



2m 

E bes t 



2m 



B ^ %l 2m 
|VS;(z)|| 2 Wz) 



Vj 3i m-—\\VS l (x)\\ 2 \Ux) 

(78) 
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We then multiply equation 1|77[) by ipi{ x ), integrate over 
the sets D?, sum the contributions of all sets and con- 
clude that 



E 2{P°best) - E > 



best 



2m 



fl&Pbest( X ) dx 



Indeed, for the case of the quartic potential, one com- 
putes E2(pl est ) — 1.372 which is seen to be lower than 
E best = 1.404 but higher than E 1 {p° best ) = 1.342. The 
relations l|69[l . (|72|l and (|79(l combined give 



e < E x {pl est ) < E 2 {p° best ) < E, 



best i 



(80) 



which proves our previous assertion that the asymptotic 
density matrix predicted by LVP has a lower energy than 
the one given by GBF. In fact, if Ebest < °°, the last 
inequality is strict except for the case when the original 
potential and the optimized trial potential are identical. 



V. AVERAGE ENERGY AT LOW 
TEMPERATURE. THE SEMI-SUM THEOREM 

The LVP approximation is intended as a technique for 
computing finite-temperature properties of a quantum 
physical system, properties that are usually of the form 



(0)p = 



f K dxp b est(x; (3)0(x; (3) 
J R dxp best (x; (3) 



Such averages can be estimated for fairly complex sys- 
tems by Monte Carlo simulations [3^. The problem we 
address in this section is the low temperature limit of 
different energy estimators. For operators which are di- 
agonal in the coordinate representation, for example the 
potential energy V(x), the estimating function 0{x) does 
not depend upon temperature and the zero temperature 
limit is 



lim (0) p 



k dx P°best( X ) ( X ) 

Ldxpl est (x) 



In this paper, we assume that the pointwise optimiza- 
tion in the configuration space involved by LVP can 
be rapidly implemented by standard local optimization 
procedures, iterative methods or other approximations. 
Since this is a big assumption by itself, estimators ex- 
plicitly depending upon the derivatives of the optimizing 
parameters B(x, x'; /?) are clearly out of question. In the 
remainder of this section, we shall consider the impor- 
tant problem of computing the ensemble average energy 
with the help of the so-called T, respectively estimat- 
ing functions, both of which can be put in a form that 
satisfies the aforementioned restriction. With regard to 
the zero temperature limit, we are interested to learn 
whether we can recover fully or only partially Ei(p best ), 
the best energy predicted by LVP. We shall prove that 
we recover Eb es t with the help of the T-method estima- 
tor, respectively the semi-sum of E\{p best ) and E 2 (pl est ) 



by the H-method estimator. The last fact is called the 
semi-sum theorem. 

We begin by considering some preliminary results. The 
maximum condition (|20ll implies the equality 



d_ 

1Tb 



p a b (x,x';l3) 



Mx.x'e 



51 



b=B(x,x';0) 



Another consequence of the same extremum condition is 
that the hessian matrix 



db 2 



p- b {x,x'-p) 



(82) 



b=B(x,x';/3) 



is negative definite Vs,x' S R. Moreover, the symme- 
try of p? (x, x'; 0) in the arguments x and x' implies the 
symmetry of the minimizing functions B(x,x';(3) in the 
same arguments. We then have the equality 



d_ 

dx 



B(x,x';f3) 



_0_ 

dx 



7 B(x,x';p) 



(83) 



At finite temperature, because of the thermal averag- 
ing, it is safe to assume that the optimizing parameters 
B (x, x'\ /3) are nice functions in their arguments with con- 
tinuous partial derivatives at least to the first order. This 
might not be true for B(x, x'; oo) which may be constant 
on the interior of some sets Di, but vary suddenly at their 
frontier. 

For the rest of the paper, we shall assume 
that p begt (x,x') is in the Sobolev space 7J 1,2 (]R 2 ). Thus, 
the norms (defined here by their square) 



\Pbest\ 



dx dx'p° best (x,x') 



'\2 



(84) 



and 



P°bestf + // [(dxP° best ) 2 + (dx'P°b est f 



\Pbest\\l 

.1 J - 

(85) 

are finite (for the case analyzed in the previous section, 
it is rather trivial to prove that these conditions are ful- 
filled). We shall also assume the existence of the sec- 
ond derivatives almost everywhere. Some mathematical 
difficulties force us to restrict the analysis to potentials 
bounded from below — thus, positive by a change of ref- 
erence. 

To avoid the excessive use of big vertical lines, we shall 
follow the rule that all functions f(x,x',b;(3) explicitly 
depending upon b are evaluated at b = B(x, x'; (3) if the 
results are to hold. The evaluation is done after any 
differentiation but before integration. 

There are a couple of energy estimators in litera- 
ture |33| , of which we shall consider the most important 
two: the so called T-method and 7i-method estimators. 
The first one is computed by temperature differentiation 
of the canonical partition function 



Pbest{x; P)dx 



(86) 
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With the help of (|81|l . one can show that 



J R pf(x;(3)dx 



(87) 



expression that is seemingly easier to compute since it 
does not involve the evaluation of the partial derivatives 
of B(x, x'; P) with respect to temperature. The low tem- 
perature limit is computed by replacing in formula 186(1 
the asymptotic density matrix given by l(51() , to produce 



lim (E) g = E best 



(88) 



i.e., the groundstate energy we get by using the T-mcthod 
estimator coincides with the best energy provided by the 
analog centroid based approximations. 

In the particular case of the HO-LVP approximation, 
the diagonal density matrix takes the form: 



2-kK 2 ^ 



hoifiC) 



x exp 
x exp | — 
where 



^p 3 B 2 h 5 (pC) + ^ 2 C 2 h 4 (pC) 



(3 I dt V t ,u 
o 



x-aP 3/2 Bh 2 (PC, t) 



2(3h 2 



(89) 



(90) 



The T-estimator function has the expression 



1 



dt V t 



x-ap 3/2 Bh 2 {(3C, t) 



-ct/3 3/2 B / dtV t „ x-a{3 3/2 Bh 2 ({3C,t) h 2 (f3C,t) 
2 Jq ' L J 



x - af3 z/2 Bh 2 {f3C,t)\hi{f3C,t) (91) 



The ii-method estimator is the direct expected value 
of the Hamiltonian 

J R dxHp best (x,x';0)\ 
J R dxp best (x; p) 

In computing the kinetic term of ((92(1 . the following for- 
mula proves beneficial: 

- H 2 k& X (d^-&)pbest{x,x'-,l3) 



(K) 



4m 



J R dxp best (x; f3) 



(93) 

We compute the expected kinetic energy as the semi-sum 
of the two identical terms, because this way no deriva- 
tives of B{x, x'; 0) appear in the final formula. By differ- 
entiation of ((81(1 against x' , we get the system of equa- 
tions 

' )2 7 pl(x, x'; P)+-^pl(x, x';(3) dB(x ^' ;(]) = (94) 



dbdx 



dx' 



and there is a similar one for the derivatives against x. 
From 1(81(1 and 1(94(1 . the following equalities can be de- 
duced by explicit calculation: 



d 2 d 2 

7 p best (x,x';P) = Q^jP- b ( x > x '^) 

dB(x,x';f3) dB(x,x';/3) 



dxdx 
d 2 



and 



dx 



d 2 



dx 1 



(95) 



—^p best (x,x';fJ) = —^p b (x,x';f3) 



dx 2 

d 2 



dx 2 

dB{x 7 x';(3) dB(x,x';/3) 



dx 



dx 



(96) 



By adding l(9"5l) and (|9"BTl and using , we get the equal- 



ity 



d 2 



dxdx 
d 2 



jp bes t(x,x'; (3) 



[ d 2 

^p t ( x ,x'-,P)- 1 ^pl(x,x';P) 



(97) 



Relation ((97() shows that there is no need for the par- 
tial derivatives of the optimizing parameters B{x,x')0) 
against x or x' in order to evaluate the ensemble average 
energy by the H-method estimator. We shall introduce 
two additional kinetic energy estimators which serve as 
intermediate tools in our computation: 



2 k Ax T&?PK x ' x 'iP) 



and 



(K)f 2 = 



2m J R dxp best {x;(3) 
H 2 k Ax &P a b( x ^ x '^) 



2m 



J M dxp best (x; /?) 



(98) 



(99) 



and denote the respective energy estimators, obtained by 
adding the ensemble average potential energy by (E)^' 1 

and (E)p' 2 respectively. The second estimator, called in 
this work of "type 2," is always greater than the first, 
which is called of "type 1." Indeed, from ((95(1 and l(96|) 
one learns that 



(K^-Wf^iE)^ 2 



(E) 
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o f Hr 92 n a (T- R\ 9B(x,x'-l3) dB(x,x';f3) 
H JVL db 2 Pb y dx dx 



c - > 0(100) 



m J R dxp best {x; j3) 

the inequality following from the negative definitiveness 
of the hessian matrix l(82|) . Now, relation 1(97(1 says that 



(K) 



H 



and so, 



(E) 



H 



l r 

2 

1 

2 



(K) 



(E)f 1 + (E) 



H,2 



(101) 
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The reader might have already realized that the for- 
mula l|102|) is the key to the semi-sum theorem announced 
at the beginning of the section. It also implies that there 
is no need for the partial derivatives of the optimizing pa- 
rameters, in order to compute the H-method estimator. 
For the case of the HO-LVP, the H-estimator function 
has the expression 

E* u (x;f3) = i + V{x) + /3 3 C%(/?C) 
2(3 



f\ 77" 

dt V ttU 

o 



x - a(3 3/2 Bh 2 {(3C, t)] h 7 {(3C, t), (103) 



where a is defined by Eq. |9*U1 

To continue, we turn our attention to some conver- 
gence problems. The expected values of the potential en- 
ergy or other diagonal operators (including the constant 
functions) converge smoothly to the expected values com- 
puted with respect to p best (x) as j3 — > oo, fact that can 
be justified in most cases with the help of the Dominated 
Convergence Theorem (see Theorem 2.24 of Ref. 134^1 . 
However, this theorem cannot be used directly in the case 
of the expected values of the operators whose estimators 
explicitly involve the partial derivatives of B(x, x'; (3). As 
we saw in the previous section, their moduli may blow 
up on dDi so no dominating function might exist. 

Using the asymptotic form 14211 predicted by Theo- 



rem 3, one easily proves the following analog of i|79|) 
h 2 f M Mx) 2 \\VS- b (x)\\ 2 dx 



(E) 



H,2 



E, 



best 



2m 



(104) 



Since the expression on the right-hand side does not ex- 
plicitly involve derivatives of the optimizing coefficients, 
we have the equality 



lim (E)f 2 



best 



^ J R P° best (x)\\VS b (x)\\ 2 dx 



2m 



(105) 

where we set b = B(x,oc) before integration. With the 
hypothesis that the potential is positive in mind, we see 
that (E)p' 2 is the biggest estimator around, so all above 
defined estimators are bounded by E best for low enough 
temperature. Moreover, the estimator (K)^' 1 cannot be 
negative since 



lim (K) 



H,l 

P 



n 2 JrIIW^IIVb^oo)^ 



2m 



lRp°best( X ) dx 



> 0. 



(106) 

It follows that the first term from the relation (|100|l is 
bounded by E best when (3 is large and again with the 
help of Theorem 3, one may establish the result 



^E best / p° best (x)dx > lim 



Mi = lim 



dx^ pt (x-,P) d ^^ dS ^ X '^ 



dx 



db 2rbK lr ' dx dx 
d 2 E{4>l) d 2 M x ) 2 1 9B(x, x';(3) dB(x, x'; 0) 



> M 2 



dx lim 



db 2 



db 2 



dx 



dx 



(107) 



d 2 E{4>D d 2 U x ? 



db 2 



db 2 



dB(x,x';f3) dB(x,x';/3) 



dx 



dx 



where we used Fatou's lemma for the last inequality. In 
this respect, notice that the evaluation of the hessian 
matrices is done at b = B(x, (3) and that 



d 2 E{4>D d 2 M x ) 



db 2 



db 2 



(108) 



is positive definite at each point x. In order for the in- 
equality to hold for arbitrarily large (3, when the energy 
hessian matrix also becomes positive definite, the follow- 
ing condition is necessary: 



d 2 E(cj>l) dB{ x, x'; oo) dB(x, x'; oo) 



db 2 



dx- 



dx 



a.e 



(109) 

or otherwise, the argument of the last integral in H107(l 
becomes arbitrarily large on a set of strictly positive 
measure, which contradicts the fact that the integral is 



bounded on the entire low temperature range. The equal- 
ity (|109|) can be realized cither by the a.e. vanishing of 
the derivatives of the optimizing coefficients as in the case 
studied in Section III, or by the vanishing of some normal 
modes of the energy hessian matrix. Consequently, 



lim c 0E^ t .d 2 dB(x,x';(3) dB(x,x';(3) 

0™ e db 2Ph[ ' P) dx dx 



d 2 ip b (x) 2 dB(x, x'; oo) dB(x, x'; < 



db 2 



dx 



dx 



>0 (110) 



and the last function integrates to M 2 < Mi. Now, we 
have 



lim 

P^oc 



,PE b . 



ff 2 



dxdx 
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pointwise, but also 
da; 



lim 



= / dx- 



dxdx 



d 2 



dxdx 

A comparison with (|95|l produces 

d 2 



(112) 



lim 



dx 



dx 

o 2 



o0E„ 



dxdx 



dxdx 

~,Pbest( x i x ) 



7 p best (x,x';(3) 

- x—x' 

+ Mi - M 2 . (113) 



With these results at hand, it is not hard to conclude 
that 



H 



1™ ( E )p 

p— ►oo K 



lim (E) 



Mi 



2m kP°best( X ) Ax 

h 2 Mi - M 2 



2m kPbest( X ) dx 

In an analog manner but using l(96|l . one proves 

h 2 Mi 



(114) 



1™ (^)^ 

p— >oo ^ 



lim (£7) o : 

P ^I K P°be S t( x ) dx 



E z(Pbest) 



Mi - M 



^2 



2m Ja^t^)^ 

Summation of 1)114(1 and (j 1 1 5(1 produces the theorem: 

Theorem 5 (Semi-sum) TTie low temperature limit of 
the H-method energy estimator is: 



lim = 1 [£i(P& esi ) + E 2 (p° best )} . (116) 

p— *OC ^ Z 



Because Mi > M 2 , the various estimators introduced 
in this section can be put in the following order: 



lim (E)f 1 <E 1 (p° best )< lim (E) H 

(j—yoo p — *oo 



< E 2 {p° bESt ) < lim {E)f/ < lim (E)^ (117) 

p— >oo ' p— >oo ^ 

[for the last inequality use (|88|l and (l|lU5p ] . For continuity 
reasons, it is convenient to define the energy of the LVP 
groundstate density matrix as 



E (Pbest) = 2 i E l(Pbest) + E 2(p°best)} ■ 



(118) 



If the decomposition in sets Di is true, the almost every- 
where vanishing of the derivatives of the optimizing coef- 
ficients implies M 2 = and we recover (|79|l as it should. 
But in this very important case, maybe more significant 
is the fact that E i{pl est ) and E (pl est ) are above the true 
groundstate energy and therefore LVP is able to provide 
a variational energy which is better than Eb est , the best 
energy predicted by the variational centroid based tech- 
niques. 



VI. THE FREE PARTICLE REFERENCE CASE 
AS THE BASIC PROTOTYPE 

To summarize the results obtained in this paper in 
one sentence, the dependence of the density matrix with 
the coordinates x and x' is reproduced by the LVP ap- 
proximation in a significantly better way than the depen- 
dence with the inverse temperature (3. This is why the 
H-method estimator behaves in a better way than the T- 
mcthod estimator. It is then interesting to compare the 
variational centroid method with the LVP method for the 
simple case when the reference system is the free particle 
one, so that V'(x) = 0. The point is that in this case we 
can leave any parameter optimization issues aside. 

Letting a = y/WfJ/m, the LVP density matrix approx- 
imation takes the form 



Po^Qr,x';/3) 
p fp (x,x';(3) 



cnp { 3 J EV x{t)+aB° t dt 
i'Np| -3 I F t , [z(i)]dtJ , (119) 



where 

v t ,o(y) 



exp 



2T 2 (t) 



V(y + z)dz, 



(H5) with r§(t) defined by 



rg(t) = a 2 E(B t ) 2 =a 2 t(l-t). 

This is the zero order approximation of the so-called par- 
tial averaging method |24| . 

The variational centroid expression for the diagonal 
centroid density matrix is Q 



p c {x-(3) = {2na 2 )- 1 ' 2 e'^ s \ 



(120) 



where 



K[y) = [ 



r y/2ira 2 /12 



exp 



2a 2 /12 



V(y + z)dz. 



The question we want to answer is which one of the 
formulae (|119|l and (|120fl provides a better description of 
the physical system. To this end, notice that the spread 
in the partial averaging formula is on average twice as 
large as the one for the centroid approximation: 



Tl(t)dt = cr 2 /6 = 2(ct 2 /12). 



This is so because the centroid position is defined as a 
path average, being the unique value x around which the 
fluctuation J (B^ — x) 2 dt of a path is minimized. It is 
therefore expected that the centroid approximation be- 
haves in a better way than the zero order partial averag- 
ing formula as far as the "direct" finite temperature par- 
tition function and the related T-method estimator are 
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concerned. This should be true even if the analysis per- 
formed in this paper showed that the high temperature 
and the low temperature limits are the same for the two 
methods. However, this does not mean that the centroid 
formula gives the better description of the system. To the 
contrary, we assert that by means of the H-method esti- 
mator, the zero order partial averaging formula provides 
the better description of the system as far as the average 
energy (and by integration against temperature, the ra- 
tio of the partition functions at different temperatures) is 
concerned. We shall present numerical evidence support- 
ing our claims by analyzing a simple case of a periodic 
monodimensional potential. We choose a periodic po- 
tential because in this case the low temperature limits of 
both the partial averaging and the centroid formulae are 
well defined. Nevertheless, the reader should be aware of 
the fact that the free particle reference is the worst sce- 
nario for LVP as to its advantage over the equivalent cen- 
troid approximation. For the HO-LVP theory, the value 
of rj(t) is controlled by the spread of the best fitting 
Gaussian and to a less extent by the temperature. Even- 
tually, for low enough temperature, L T^(t)dt equals the 
spread of the best fitting Gaussian, but this also happens 
for the centroid based approximations. Therefore, the 
latter's advantage is diminished. 

Let us consider a monodimensional periodic potential 
of period 2L and let 



V(x) = 



v k e 



ikTTX j L 



(121) 



be its Fourier series. By the reality of the potential V(x) 



we have v- 



By Theorem 3, the low temperature 



asymptotic of the zero order partial averaging density 
matrix is 

p£ A (x, x'; 13) a exp{-[S(:r) + S(x')}} exp(-0vo)/>/2L, 

(122) 

where 



2L 



V(x)dx 



is the cell average of the potential and where 
2mL 2 



S(x) 



^fe ikTrx j L 

h 2 ir 2 ^ k 2 



E 



(123) 



Then, the low-temperature limit of the T-method esti- 
mator for both the centroid and the partial averaging 
approximations is 



lim (E) T „ = v , 

(3— »oo ^ 



(124) 



while the low temperature limit of the H-method estima- 
tor for the partial averaging formula is 



lim (E) 



^0 



h 2 J_ L eyvs>[-2S{x)\ \\\7S(x)\\ 2 &x 



2m 



f_ L exp[-2S(x)]da 



(125) 
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FIG. 4: The groundstate energies given by the low temper- 
ature limits of the T-method estimator and the H-method 
estimator are plotted together with the exact groundstate en- 
ergies for various values of the potential frequency v. 



For a periodic potential, a state is called bounding if 
its expected energy is strictly smaller than the potential 
average. If follows that the centroid method fails to pre- 
dict a bounding groundstate for the periodic potential. 
To the contrary, the zero order partial averaging formula 
does predict a bounding groundstate by means of the H- 
estimator. The predicted groundstate energy is shown in 
Fig. 4 for the case of the periodic potential 

V(x) = 0.5 [1 + cos(2twx)] 

and is plotted as a function of the frequency v. Again we 
use atomic units and consider a particle of mass m = 1 . 
The exact groundstate energies were computed with the 
help of the Rayleigh-Ritz variational principle by expan- 
sion in a Fourier series. One notices that the H-method 
energy is in good agreement with the exact result in the 
range of high frequencies, but extrapolates to vq/2 in- 
stead of in the low frequency range. 

It is instructive to study the behavior of the finite tem- 
perature energy estimators for the periodic potential (we 
set v = 0.25). Because in practice one stops the calcula- 
tions the first time the energy fails to decrease with the 
temperature, Fig. 5 suggests that the centroid T-method 
estimator has a better behavior than the corresponding 
LVP T-method energy, even if their low temperature lim- 
its are the same. However, the LVP H-method estimator 
has an overall even better behavior, providing the closest 
answer to the exact energy for the investigated range of 
temperatures. The exact energies were computed by nu- 
merical matrix multiplication (see for instance, Ref. |35j). 



VII. CONCLUDING REMARKS 

The Local Variational Principle is alleviating some of 
the disadvantages of the GBF principle. Specifically, it 
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Values of the partial averaging T-estimator Ea , par- 

and 



tial averaging H-estimator Ep, centroid T-estimator Ep 
the exact energy E@ x as functions of the inverse temperature 
P in a logarithmic scale. Atomic units were used for energy 
and p. 



makes better use of simple reference potentials to pro- 
duce accurate diagonal density matrices. Moreover, the 
Trotter composition method systematically improves the 
finite temperature LVP density matrix up to the cor- 
rect value. Then, LVP always gives a better groundstate 
energy than the one predicted by the Gibbs-Bogoliubov- 
Feynman principle provided that the H-method estima- 
tor is used. Finally, we conclude that LVP gives a de- 
scription of the physical system which is more accurate 
and more complete than the one provided by the centroid 
based approximations. I remind the reader that this sit- 
uation was caused by the fact that there is no practical 
way of defining an H-estimator for the latter methods. 
Since the GBF principle becomes essentially a local ap- 
proximation in the centroid space, it may be possible that 
by modifying the centroid idea, one might be able to con- 
struct a new technique inheriting the best features of the 
two methods. 

Unfortunately, there are several drawbacks of the LVP 
method (all of them being shared by the variational cen- 
troid techniques). The first problem is the need for ex- 
pensive pointwisc local maximization procedures. How- 
ever, they can be approximately performed as long as 
the estimators of the evaluated properties do not explic- 
itly involve the derivatives of the optimizing parameters 
(both the T-method and the H-method estimators are 
"stable" in this respect). To understand the second unde- 
sirable feature, let us imagine that we slightly unbalance 
the double well potential. Then, the best Gaussian dis- 
tribution that realizes the minimum of l|48l) is unique and 
will be localized in the well of lower energy. Thus, there 
is a discontinuity of the method with respect to the bal- 
ancing of the main wells of the potential. LVP partially 
accounts for this problem with the help of a correction 
factor (see Eg. I44f) . but the improvement is not always 
satisfactory. 



Finally, let us notice that given a point (x,x') in the 
configuration space, the HO-LVP technique gives explicit 
information about how the density matrix looks around 
this point on a range established by the temperature and 
by the quantum properties of the system. This informa- 
tion can be used to create optimal filters in conjunction 
with standard semiclassical approximations for the ther- 
malized quantum dynamical correlation functions. In 
this respect, as discussed by Miller [3(|, the semiclassi- 
cal Van Vleck or Herman-Kluk propagators in the initial 
value representation could be effectively used in quantum 
simulations for sufficiently large systems and for chemi- 
cally relevant times if not for the associated "sign prob- 
lem." This is generated by the integration over the initial 
conditions and can be alleviated by use of filtering tech- 
niques, which are also advantageous because they require 
only local information about the density matrix for an 
optimal implementation. Though a final decision as to 
the feasibility of this idea awaits future detailed work, I 
anticipate that this local (analytical) information can be 
quantitatively furnished by the HO-LVP technique under 
the form of a "built in" filter. 
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APPENDIX A: DEFINITION OF THE 
H-FUNCTIONS 
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